Ejercicio 1¶

1. Reporte descriptivo del dataset¶

Dimensiones del dataset, número de variables continuas y categóricas, distribuciones y comentarios generales

In [1]:
import pandas as pd
import numpy as np

pd.set_option('display.max_rows', None, 'display.max_columns', None)

fev = pd.read_csv(r'FEV_data.csv')

fev.head()
Out[1]:
seqnbr subjid age fev height sex smoke
0 1 301 9 1.708 57.0 2 2
1 2 451 8 1.724 67.5 2 2
2 3 501 7 1.720 54.5 2 2
3 4 642 9 1.558 53.0 1 2
4 5 901 9 1.895 57.0 1 2
In [2]:
fev.info()
fev.shape

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 654 entries, 0 to 653
Data columns (total 7 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   seqnbr  654 non-null    int64  
 1   subjid  654 non-null    int64  
 2   age     654 non-null    int64  
 3   fev     654 non-null    float64
 4   height  654 non-null    float64
 5   sex     654 non-null    int64  
 6   smoke   654 non-null    int64  
dtypes: float64(2), int64(5)
memory usage: 35.9 KB
Out[2]:
(654, 7)

Después de cargar el dataset y ver la información del mismo, vemos que tenemos 654 observaciones y 7 columnas, sin embargo la primera (seqnbr) es unicamente el número de secuencia o index de cada observación, mientras que subjid es el ID del paciente (en la buena teoría debería ser un número random). Por otro lado las columnas de sex y smoke son categóricas (realmente binarias), por lo que vamos a convertirlas al tipo 'category', no sin antes asegurarnos de que efectivamente contienen unicamente dos valores diferentes cada una.

In [3]:
fev.nunique()
Out[3]:
seqnbr    654
subjid    654
age        17
fev       575
height     56
sex         2
smoke       2
dtype: int64
In [4]:
factors = list(fev.loc[:, fev.nunique()<= 10])
fev[factors] = fev[factors].astype('category')
In [5]:
fev.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 654 entries, 0 to 653
Data columns (total 7 columns):
 #   Column  Non-Null Count  Dtype   
---  ------  --------------  -----   
 0   seqnbr  654 non-null    int64   
 1   subjid  654 non-null    int64   
 2   age     654 non-null    int64   
 3   fev     654 non-null    float64 
 4   height  654 non-null    float64 
 5   sex     654 non-null    category
 6   smoke   654 non-null    category
dtypes: category(2), float64(2), int64(3)
memory usage: 27.2 KB

Con esto, tendríamos una variable objetivo (fev) y 4 son posibles variables predictoras, de las cuales tenemos 2 variables continuas (age y heigth) y 2 variables categóricas.

No tenemos missing values puesto que todas las columnas tienen 654 datos, sin embargo debemos comprobar que no haya outliers, además queremos visualizar la distribución de las variables. Para esto haremos una tabla con datos descriptivos y posteriormente una inspección gráfica.

In [6]:
fev.describe()
Out[6]:
seqnbr subjid age fev height
count 654.00000 654.000000 654.000000 654.000000 654.000000
mean 327.50000 37169.571865 9.931193 2.636780 61.143578
std 188.93782 23690.860350 2.953935 0.867059 5.703513
min 1.00000 201.000000 3.000000 0.791000 46.000000
25% 164.25000 15811.000000 8.000000 1.981000 57.000000
50% 327.50000 36071.000000 10.000000 2.547500 61.500000
75% 490.75000 53638.500000 12.000000 3.118500 65.500000
max 654.00000 90001.000000 19.000000 5.793000 74.000000
In [7]:
fev.describe(exclude = np.number)
Out[7]:
sex smoke
count 654 654
unique 2 2
top 1 2
freq 336 589
In [8]:
import seaborn as sns
import matplotlib.pyplot as plt
import plotly.express as px

def histogram_boxplot(data, xlabel = None, title = None, font_scale=2, figsize=(9,8), bins = None):
    # Crear ventana para los subgráficos
    f2, (ax_box2, ax_hist2) = plt.subplots(2, sharex=True, gridspec_kw={"height_ratios": (.15, .85)}, figsize=figsize)
    # Crear boxplot
    sns.boxplot(x=data, ax=ax_box2)
    # Crear histograma
    sns.histplot(x=data, ax=ax_hist2, bins=bins) if bins else sns.histplot(x=data, ax=ax_hist2)
    # Pintar una línea con la media
    ax_hist2.axvline(np.mean(data),color='g',linestyle='-')
    # Pintar una línea con la mediana
    ax_hist2.axvline(np.median(data),color='y',linestyle='--')
    # Asignar título y nombre de eje si tal
    if xlabel: ax_hist2.set(xlabel=xlabel)
    if title: ax_box2.set(title=title, xlabel="")
    # Mostrar gráfico
    plt.show()


def plot(col):
     if col.dtypes != 'category':
        print('Variable Continua')
        histogram_boxplot(col, xlabel = col.name, title = 'Distibución continua')
     else:
        print('Cat')
        fig = px.bar(col.value_counts())
        fig.show()
        


fev.drop(['seqnbr', 'subjid'], axis = 1).apply(plot)

Variable Continua

Variable Continua

Variable Continua

Cat

Cat
Out[8]:
age       None
fev       None
height    None
sex       None
smoke     None
dtype: object

No tenemos missing data, no vemos errores en los datos; parece haber algunos outliers en la variable fev, sin embargo al tratarse de la variable objetivo y además, al no estar tampoco tan alejados de los demás valores al punto de parecer irreal, seguiremos trabajando con ellos para no perder posible información valiosa.

En cuanto a la distribución, tenemos para age y height distribuciones normales, fev, no es tan simétrica, pero de igual forma, posiblemente esté muy cerca de la normalidad. Mientas que las variables categóricas, tenemos por un lado Sex, que está muy bien distribuida entre sus dos categorías, con cantidad de datos muy similar para cada una; y por el otro lado, smoke (que como sería de esperarse al tratarse de niños y adolescentes) tiene muy poca representación de pacientes fumandores, la mayoría está en la categoría 2 (asumiremos que es no fumadores).

Parece ser una base de datos bastante limpia y lista para modelar!

2. Decidir si se descarta de inicio alguna de las variables de cara al modelado¶

Como se mencionó antes, la variable seqnbr, es el número de secuencia, por lo que no trabajaremos con ella.

La variable subjid pareciera ser también algun número de identificación del paciente, en ese caso no debería tomarse en cuenta puesto que se trata de números aleatorios.

La variable smoke, está muy poco representada en una de sus dos categorías, por lo que podría aportarnos muy poca información e incluso información sesgada, no vamos a descartarla desde ya, sino que trabajaremos con ella para ver si de alguna forma nos aporta información valiosa, pero está en la mira para descartarla en caso de que no nos aporte mucho.

3. Modelo de Regresión lineal sin interacciones entre las variables¶

Vamos a hacer una inspección visual y general, para hacernos una idea de las posibles relaciones entre las variables continuas

In [9]:
g = sns.PairGrid(fev.iloc[:, 2:])
g.map_upper(sns.scatterplot)
g.map_diag(sns.histplot)
g.map_lower(sns.regplot)
Out[9]:
<seaborn.axisgrid.PairGrid at 0x1c818948250>

Parece haber una clara relación entre las variables de edad y altura con la variable objetivo y también parece haber una dependencia entre ellas; sin embargo aún no hemos visto nada con respecto a sex y smoke que son categóricas.

Utilicemos pandas profiling para tener más información:

In [10]:
import pandas_profiling
pandas_profiling.ProfileReport(fev.iloc[:, 2:])

Summarize dataset:   0%|          | 0/5 [00:00<?, ?it/s]
Generate report structure:   0%|          | 0/1 [00:00<?, ?it/s]
Render HTML:   0%|          | 0/1 [00:00<?, ?it/s]
Out[10]:




























En este reporte vemos resumido mucho de lo que analizamos antes, sin embargo hay datos adicionales en la sección de correlaciones: por un lado nos confirma una correlación fuerte entre age y heigth con la variable objetivo (y entre ellas también!) y por otro lado, nos indica que hay cierta relación entre sex y smoke con fev, aunque no es tan fuerte. Para esta segundo punto (correlación categórica-numérica) el reporte utiliza el método v de cramer; por lo que sería interesante ver este análisis en detalle.


Primero vamos a separar la variable objetivo para poder manipular las demás variables sin riesgo de modificar o transformar la variable de respuesta



















In [11]:



     


varObj = fev.fev
inputs = fev.drop(['seqnbr', 'fev', 'subjid'],axis=1)
inputs.info()



     






















    
    






<class 'pandas.core.frame.DataFrame'>
RangeIndex: 654 entries, 0 to 653
Data columns (total 4 columns):
 #   Column  Non-Null Count  Dtype   
---  ------  --------------  -----   
 0   age     654 non-null    int64   
 1   height  654 non-null    float64 
 2   sex     654 non-null    category
 3   smoke   654 non-null    category
dtypes: category(2), float64(1), int64(1)
memory usage: 11.9 KB




























Y ahora si, usaremos la función de VCramer para un análisis preliminar de la relación de las variables con la variable objetivo (aunque ya tenemos una ligera idea de esto con base en la inspección visual que hicimos, pero queremos hacer la comparación de la correlación para cada input)



















In [12]:



     


import scipy.stats as stats

def cramers_v(var1, vObj):
    
    if not var1.dtypes == 'category':
        bins = min(5,var1.value_counts().count())
        var1 = pd.cut(var1, bins = 5)
    if not vObj.dtypes == 'category':
        bins = min(5,vObj.value_counts().count())
        vObj = pd.cut(vObj, bins = 5)
        
    data = pd.crosstab(var1, vObj).values
    vCramer = stats.contingency.association(data, method = 'cramer')
    return vCramer


# Ejemplo uso univariante
#cramers_v(vinosCompra['Etiqueta'],vinosCompra['Beneficio'])

# Aplicar la función al input completo contra la objetivo
tablaCramer = pd.DataFrame(inputs.apply(lambda x: cramers_v(x,varObj)),columns=['VCramer'])

# Obtener el gráfico de importancia de las variables frente a la objetivo continua según vcramer
import plotly.express as px
px.bar(tablaCramer,x=tablaCramer.VCramer,title='Relaciones respecto al FEV').update_yaxes(categoryorder="total ascending")



     






















    
    







                            
                    




























Ahora si, podemos ver ordenados los inputs según su correlación con la variable objetivo!


A lo que vinimos, vamos a generar el modelo de regresión lineal para la predicción del FEV; lo primero es separar los datos en dos grupos: training y test; posteriormente generamos el modelo con los datos de training



















In [13]:



     


from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(inputs, varObj, test_size=0.2, random_state=42)

# Comprobamos dimensiones
print('Training dataset shape:', X_train.shape, y_train.shape)
print('Testing dataset shape:', X_test.shape, y_test.shape)



     






















    
    






Training dataset shape: (523, 4) (523,)
Testing dataset shape: (131, 4) (131,)
























In [14]:



     


# Unimos los datos de training de la variable objetivo, con los datos de training de los inputs:

fev_training = X_train.join(y_train)
fev_training.head()



     






















    
    
Out[14]:











  
    

      
      age
      height
      sex
      smoke
      fev
    

  
  
    

      321
      12
      71.0
      1
      2
      4.550
    

    

      456
      12
      68.0
      1
      2
      4.411
    

    

      340
      10
      62.0
      1
      2
      1.937
    

    

      29
      9
      60.0
      2
      2
      2.100
    

    

      570
      10
      65.0
      1
      2
      3.090
    

  



























In [15]:



     


# Generamos el modelo completo:

from statsmodels.formula.api import ols 

reg_lineal = ols('fev ~ age + height + sex + smoke', data=fev_training).fit()
reg_lineal.summary()



     






















    
    
Out[15]:








OLS Regression Results


  Dep. Variable:           fev         R-squared:             0.766 




  Model:                   OLS         Adj. R-squared:        0.764 




  Method:             Least Squares    F-statistic:           423.4 




  Date:             Sun, 05 Mar 2023   Prob (F-statistic): 1.09e-161




  Time:                 09:45:19       Log-Likelihood:      -281.13 




  No. Observations:         523        AIC:                   572.3 




  Df Residuals:             518        BIC:                   593.6 




  Df Model:                   4                                     




  Covariance Type:      nonrobust                                   








                coef     std err      t      P>|t|  [0.025    0.975]  




  Intercept     -4.3576     0.271   -16.067  0.000    -4.890    -3.825




  sex[T.2]      -0.1511     0.038    -4.028  0.000    -0.225    -0.077




  smoke[T.2]     0.0652     0.066     0.983  0.326    -0.065     0.196




  age            0.0632     0.011     5.678  0.000     0.041     0.085




  height         0.1041     0.005    19.150  0.000     0.093     0.115








  Omnibus:       11.617   Durbin-Watson:         2.092




  Prob(Omnibus):  0.003   Jarque-Bera (JB):     17.038




  Skew:           0.178   Prob(JB):           0.000200




  Kurtosis:       3.810   Cond. No.               932.





Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.



























Como habíamos mencionado antes, la variable smoke está muy poco representada en una de sus categorías (es normal, los adolescentes no deberían fumar), además vemos que el p-value de su coeficiente de correlación es mayor a 0.05 (es el único mayor a 0); por lo que, como habíamos indicado antes, vamos a excluirla del modelo y ver cual es el efecto en el R cuadrado:



















In [16]:



     


reg_lineal1 = ols('fev ~ age + height + sex', data=fev_training).fit()
reg_lineal1.summary()



     






















    
    
Out[16]:








OLS Regression Results


  Dep. Variable:           fev         R-squared:             0.765 




  Model:                   OLS         Adj. R-squared:        0.764 




  Method:             Least Squares    F-statistic:           564.2 




  Date:             Sun, 05 Mar 2023   Prob (F-statistic): 6.86e-163




  Time:                 09:45:52       Log-Likelihood:      -281.61 




  No. Observations:         523        AIC:                   571.2 




  Df Residuals:             519        BIC:                   588.3 




  Df Model:                   3                                     




  Covariance Type:      nonrobust                                   








               coef     std err      t      P>|t|  [0.025    0.975]  




  Intercept    -4.2829     0.260   -16.450  0.000    -4.794    -3.771




  sex[T.2]     -0.1540     0.037    -4.117  0.000    -0.227    -0.081




  age           0.0601     0.011     5.629  0.000     0.039     0.081




  height        0.1043     0.005    19.223  0.000     0.094     0.115








  Omnibus:       11.573   Durbin-Watson:         2.100




  Prob(Omnibus):  0.003   Jarque-Bera (JB):     17.840




  Skew:           0.155   Prob(JB):           0.000134




  Kurtosis:       3.850   Cond. No.               892.





Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.



























El R cuadrado del modelo practicamente no cambió al excluir smoke, por lo que lo mantendremos así. Parece ser que tenemos un buen modelo con un R2 de 0.765; sin embargo solo por entenderlo un poco mejor, vamos a ver cual es el aporte de cada variable al R2:



















In [17]:



     


import patsy

# Generamos las matrices de diseño según la fórmula de modelo completo
y, X = patsy.dmatrices('fev ~ age + height + sex + smoke', fev_training, return_type='dataframe')



     



















In [18]:



     


from relativeImp import relativeImp

names=X.columns.tolist()[1:]

# Calculamos importancia relativa
df_results = relativeImp(X.join(y), outcomeName = 'fev', driverNames = names)

# Ordenamos valores 
df_results.sort_values(by='normRelaImpt', ascending=False)



     






















    
    
Out[18]:











  
    

      
      driver
      rawRelaImpt
      normRelaImpt
    

  
  
    

      3
      height
      0.443669
      57.937958
    

    

      2
      age
      0.277787
      36.275651
    

    

      0
      sex[T.2]
      0.024235
      3.164859
    

    

      1
      smoke[T.2]
      0.020075
      2.621533
    

  



























In [19]:



     


px.bar(df_results,x='normRelaImpt',y='driver',title='Importancia relativa por aportación al R2').update_yaxes(categoryorder="total ascending").show()



     






















    
    







                            
                    




























Según esto, podríamos eliminar (además de smoke) también a sex del modelo, veamos como cambia:



















In [20]:



     


reg_lineal2 = ols('fev ~ age + height', data=fev_training).fit()
reg_lineal2.summary()



     






















    
    
Out[20]:








OLS Regression Results


  Dep. Variable:           fev         R-squared:             0.758 




  Model:                   OLS         Adj. R-squared:        0.757 




  Method:             Least Squares    F-statistic:           812.9 




  Date:             Sun, 05 Mar 2023   Prob (F-statistic): 8.87e-161




  Time:                 09:46:25       Log-Likelihood:      -290.02 




  No. Observations:         523        AIC:                   586.0 




  Df Residuals:             520        BIC:                   598.8 




  Df Model:                   2                                     




  Covariance Type:      nonrobust                                   








               coef     std err      t      P>|t|  [0.025    0.975]  




  Intercept    -4.5960     0.253   -18.181  0.000    -5.093    -4.099




  age           0.0523     0.011     4.900  0.000     0.031     0.073




  height        0.1095     0.005    20.417  0.000     0.099     0.120








  Omnibus:       18.099   Durbin-Watson:         2.110




  Prob(Omnibus):  0.000   Jarque-Bera (JB):     31.130




  Skew:           0.231   Prob(JB):           1.74e-07




  Kurtosis:       4.103   Cond. No.               852.





Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.



























El R2 se disminuye ligeramente, pero nada considerable, por lo que tenemos dos posibles modelos lineales: uno que incluye edad, altura y sexo, y otro únicamente con edad y altura (que desde hace rato vimos que eran las que parecían tener mayor correlación).


Vamos a evaluar ambos modelos con los datos de test:



















In [21]:



     


from sklearn.metrics import mean_squared_error, r2_score

# Predicciones para test modelo completo
fev_pred1 = reg_lineal1.predict(X_test)

# Cálculo de performance
print("Mean squared error: %.2f" % np.sqrt(mean_squared_error(y_test, fev_pred1)))
# Coeficiente de determinación:
print("Coefficient of determination: %.2f" % r2_score(y_test, fev_pred1))



     






















    
    






Mean squared error: 0.40
Coefficient of determination: 0.80
























In [22]:



     


# Predicciones para test modelo completo
fev_pred2 = reg_lineal2.predict(X_test)

# Cálculo de performance
print("Mean squared error: %.2f" % np.sqrt(mean_squared_error(y_test, fev_pred2)))
# Coeficiente de determinación:
print("Coefficient of determination: %.2f" % r2_score(y_test, fev_pred2))



     






















    
    






Mean squared error: 0.41
Coefficient of determination: 0.79




























Tenemos datos muy similares para ambos modelos; los dos son aceptables, pero una variable menos nos simplificaría mucho el modelo. Veamos las interacciones:


4. Ajusta el mejor modelo de regresión lineal con interacciones entre las variables que resulten relevantes.¶






















Ya habíamos visto en nuestro análisis gráfico, que hay una aparente interacción entre age y height, que además como vimos antes, vienen siendo las variables más relevantes para la predicción del valor de FEV, por lo que vamos a tomar en cuenta esa interacción:



















In [23]:



     


reg_lin_interac2 = ols('fev ~ age + height + age*height', data= fev_training).fit()
reg_lin_interac2.summary()



     






















    
    
Out[23]:








OLS Regression Results


  Dep. Variable:           fev         R-squared:             0.783 




  Model:                   OLS         Adj. R-squared:        0.781 




  Method:             Least Squares    F-statistic:           622.7 




  Date:             Sun, 05 Mar 2023   Prob (F-statistic): 1.72e-171




  Time:                 09:48:36       Log-Likelihood:      -261.65 




  No. Observations:         523        AIC:                   531.3 




  Df Residuals:             519        BIC:                   548.3 




  Df Model:                   3                                     




  Covariance Type:      nonrobust                                   








                coef     std err      t      P>|t|  [0.025    0.975]  




  Intercept     -0.4138     0.593    -0.698  0.485    -1.579     0.751




  age           -0.4306     0.063    -6.789  0.000    -0.555    -0.306




  height         0.0410     0.010     4.013  0.000     0.021     0.061




  age:height     0.0077     0.001     7.712  0.000     0.006     0.010








  Omnibus:       17.104   Durbin-Watson:         2.071




  Prob(Omnibus):  0.000   Jarque-Bera (JB):     37.858




  Skew:          -0.053   Prob(JB):           6.02e-09




  Kurtosis:       4.314   Cond. No.           2.25e+04





Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 2.25e+04. This might indicate that there are
strong multicollinearity or other numerical problems.























In [24]:



     


reg_lin_interac1 = ols('fev ~ age + height + sex + age*height', data= fev_training).fit()
reg_lin_interac1.summary()



     






















    
    
Out[24]:








OLS Regression Results


  Dep. Variable:           fev         R-squared:             0.785 




  Model:                   OLS         Adj. R-squared:        0.784 




  Method:             Least Squares    F-statistic:           474.1 




  Date:             Sun, 05 Mar 2023   Prob (F-statistic): 1.49e-171




  Time:                 09:48:40       Log-Likelihood:      -258.17 




  No. Observations:         523        AIC:                   526.3 




  Df Residuals:             518        BIC:                   547.6 




  Df Model:                   4                                     




  Covariance Type:      nonrobust                                   








                coef     std err      t      P>|t|  [0.025    0.975]  




  Intercept     -0.5429     0.592    -0.918  0.359    -1.705     0.619




  sex[T.2]      -0.0968     0.037    -2.635  0.009    -0.169    -0.025




  age           -0.3881     0.065    -5.961  0.000    -0.516    -0.260




  height         0.0431     0.010     4.229  0.000     0.023     0.063




  age:height     0.0071     0.001     6.971  0.000     0.005     0.009








  Omnibus:       14.277   Durbin-Watson:         2.066




  Prob(Omnibus):  0.001   Jarque-Bera (JB):     28.474




  Skew:          -0.054   Prob(JB):           6.56e-07




  Kurtosis:       4.138   Cond. No.           2.26e+04





Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 2.26e+04. This might indicate that there are
strong multicollinearity or other numerical problems.



























De nuevo vemos que la presencia o no de sex en el modelo, no afecta mucho, de hecho afecta aún menos en el modelo que incluye la interacción de edad y altura (que por cierto, parece ser un modelo levemente mejor que el que no la incluye si nos basamos en el R2.


Vamos a ver la evaluación de estos modelos al comparar las predicciones con los datos de test y ya luego veremos la validación cruzada de ambos modelos:



















In [25]:



     


# Predicciones para test modelo completo
fev_pred3 = reg_lin_interac1.predict(X_test)

# Cálculo de performance
print("Mean squared error: %.2f" % np.sqrt(mean_squared_error(y_test, fev_pred3)))
# Coeficiente de determinación:
print("Coefficient of determination: %.2f" % r2_score(y_test, fev_pred3))



     






















    
    






Mean squared error: 0.39
Coefficient of determination: 0.82
























In [26]:



     


# Predicciones para test modelo completo
fev_pred4 = reg_lin_interac2.predict(X_test)

# Cálculo de performance
print("Mean squared error: %.2f" % np.sqrt(mean_squared_error(y_test, fev_pred4)))
# Coeficiente de determinación:
print("Coefficient of determination: %.2f" % r2_score(y_test, fev_pred4))



     






















    
    






Mean squared error: 0.39
Coefficient of determination: 0.81




























Efectivamente vemos que el aporte de sex al modelo, es mínimo, podemos excluirlo para simplificar el modelo sin perder información importante!























Comparación de ambos modelos por validación cruzada repetida. ¿Cuál de ellos tiene mejor comportamiento en generalización?¶


















In [27]:



     


from sklearn.model_selection import cross_val_score
from sklearn.model_selection import RepeatedKFold
from sklearn.linear_model import LinearRegression

def cross_val_lin(formula, data, seed=12345):
    y, X = patsy.dmatrices(formula, data, return_type='dataframe')
    model = LinearRegression()
    # Establecemos esquema de validación fijando random_state (reproducibilidad)
    cv = RepeatedKFold(n_splits=5, n_repeats=20, random_state=seed)
  
    # Obtenemos los resultados de R2 para cada partición tr-tst
    scores = cross_val_score(model, X, y, cv=cv)
  
    # Sesgo y varianza
    print('Modelo: ' + formula)
    print('Coeficiente de determinación R2: %.3f (%.3f)' % (np.mean(scores), np.std(scores)))
      
    #sns.violinplot(y=scores,palette='viridis')
      
    return(scores)



     



















In [28]:



     


form1 = 'fev ~ age + height'
form2 = 'fev ~ age + height + age*height'



# Creamos lista de fórmulas   
list_form = [form1,form2]
list_form

# Aplicamos a toda la lista la función creada (devuelve un dataframe pero está transpuesto)
list_res = pd.DataFrame(map(lambda x: cross_val_lin(x,fev, seed=2022),list_form))
list_res



     






















    
    






Modelo: fev ~ age + height
Coeficiente de determinación R2: 0.760 (0.030)
Modelo: fev ~ age + height + age*height
Coeficiente de determinación R2: 0.783 (0.033)










    
    
Out[28]:











  
    

      
      0
      1
      2
      3
      4
      5
      6
      7
      8
      9
      10
      11
      12
      13
      14
      15
      16
      17
      18
      19
      20
      21
      22
      23
      24
      25
      26
      27
      28
      29
      30
      31
      32
      33
      34
      35
      36
      37
      38
      39
      40
      41
      42
      43
      44
      45
      46
      47
      48
      49
      50
      51
      52
      53
      54
      55
      56
      57
      58
      59
      60
      61
      62
      63
      64
      65
      66
      67
      68
      69
      70
      71
      72
      73
      74
      75
      76
      77
      78
      79
      80
      81
      82
      83
      84
      85
      86
      87
      88
      89
      90
      91
      92
      93
      94
      95
      96
      97
      98
      99
    

  
  
    

      0
      0.787667
      0.716968
      0.799950
      0.729255
      0.761944
      0.791089
      0.776350
      0.688730
      0.788862
      0.765495
      0.772354
      0.786487
      0.769036
      0.718342
      0.762030
      0.784979
      0.741495
      0.785038
      0.709679
      0.772243
      0.750218
      0.789058
      0.750301
      0.784801
      0.719166
      0.711838
      0.761917
      0.725076
      0.818966
      0.767828
      0.723102
      0.712695
      0.774206
      0.809178
      0.785401
      0.823126
      0.691026
      0.811060
      0.708003
      0.747786
      0.797421
      0.736606
      0.727761
      0.768323
      0.780752
      0.777619
      0.757366
      0.727513
      0.764900
      0.772503
      0.802393
      0.762426
      0.778155
      0.698572
      0.755787
      0.803356
      0.710792
      0.782174
      0.763558
      0.727841
      0.809943
      0.768817
      0.756976
      0.748240
      0.746915
      0.73496
      0.757857
      0.791671
      0.743574
      0.768408
      0.794476
      0.765307
      0.778449
      0.736414
      0.724770
      0.760133
      0.738737
      0.783343
      0.806688
      0.720035
      0.726006
      0.721742
      0.775253
      0.781166
      0.775853
      0.731595
      0.778746
      0.72327
      0.784665
      0.798909
      0.771248
      0.750733
      0.771886
      0.727637
      0.769702
      0.775609
      0.761466
      0.763527
      0.750511
      0.760841
    

    

      1
      0.808889
      0.743159
      0.815954
      0.749891
      0.798461
      0.811031
      0.812616
      0.721601
      0.811489
      0.762009
      0.801723
      0.817343
      0.795215
      0.740051
      0.777243
      0.813767
      0.770766
      0.816728
      0.697378
      0.800250
      0.794001
      0.801966
      0.779349
      0.807909
      0.730690
      0.716065
      0.801887
      0.733507
      0.845703
      0.796791
      0.737427
      0.746988
      0.790016
      0.819267
      0.825149
      0.840542
      0.706053
      0.840481
      0.726434
      0.780371
      0.834036
      0.756675
      0.741305
      0.791445
      0.800239
      0.813321
      0.780701
      0.724106
      0.796349
      0.794709
      0.810594
      0.769931
      0.805100
      0.734141
      0.790558
      0.836799
      0.744945
      0.803868
      0.771287
      0.752825
      0.822430
      0.779777
      0.786041
      0.788226
      0.765887
      0.75224
      0.770672
      0.815152
      0.788136
      0.789405
      0.821415
      0.789580
      0.812760
      0.758850
      0.733953
      0.779442
      0.765022
      0.805122
      0.824413
      0.753467
      0.731746
      0.748743
      0.796179
      0.817650
      0.799745
      0.743433
      0.785447
      0.77693
      0.823026
      0.790566
      0.796992
      0.751288
      0.806779
      0.759025
      0.791933
      0.797006
      0.793718
      0.756039
      0.793153
      0.785192
    

  



























In [29]:



     


results = list_res.T.melt()
results.columns = ['Modelo','R2']
results.groupby(['Modelo']).describe().T



     






















    
    
Out[29]:











  
    

      
      Modelo
      0
      1
    

  
  
    

      R2
      count
      100.000000
      100.000000
    

    

      mean
      0.760046
      0.782857
    

    

      std
      0.030248
      0.032920
    

    

      min
      0.688730
      0.697378
    

    

      25%
      0.736050
      0.756516
    

    

      50%
      0.764229
      0.790287
    

    

      75%
      0.781418
      0.807061
    

    

      max
      0.823126
      0.845703
    

  































Vemos un mejor coeficiente de determinación para el modelo con interacciones, aunque también una leve desmejora en el error estandar, pero esta desmejora es mínima, entonces nos vamos a dejar influenciar más por el coeficiente de determinación (que tampoco tiene diferencia abismal, pero si se ve mejoría). Veamos de una forma gráfica la diferencia en el R2:



















In [30]:



     


sns.boxplot(x='Modelo',y='R2',data=results,palette='viridis')



     






















    
    
Out[30]:








<AxesSubplot:xlabel='Modelo', ylabel='R2'>










    
    



































Por lo que podemos ver, ambos modelos son relativamaente similares, aunque el modelo con form2 (con interacción) es el de mejor ajuste por lo que este sería nuestro modelo final!


form2 = 'fev ~ age + height + age*height'